# Alex Gazmararian
# agazmararian@gmail.com

library(tidyverse)
library(tidylog)
library(readxl)
library(here)
library(modelsummary)
library(sf)
library(tigris)
library(terra)

# Output file path
output_file <- here("data", "inter", "highway_counties_processed.csv")

# Cache the RAW tigris download, not the processed output
cache_file <- here("data", "cache", "tigris_counties_2018.rds")

# Step 1: Get county boundaries (from cache or download)
if (file.exists(cache_file)) {
  message("[OK] Using cached tigris counties from data/cache/")
  counties.pre <- readRDS(cache_file)
} else {
  message("Downloading county boundaries from tigris...")
  counties.pre <- tigris::counties(cb = TRUE, year = 2018)
  
  # Cache the raw download
  dir.create(dirname(cache_file), recursive = TRUE, showWarnings = FALSE)
  saveRDS(counties.pre, cache_file)
  message("[OK] Cached tigris counties to: ", cache_file)
}

# Step 2: Process (always runs, uses cached or fresh download)
message("Processing highway intersection data...")

# Transportation access: proximity to major highways/railroads
roads <- st_read(here("data", "input", "primary_roads", "tl_2023_us_primaryroads.shp"))

counties.pre <- st_transform(counties.pre, crs = st_crs(roads))
names(counties.pre) <- tolower(names(counties.pre))
counties.pre <- counties.pre %>%
  mutate(fips.pre = paste0(statefp, countyfp)) %>%
  dplyr::select(fips.pre)

# Subset to only interstate roads
roads <- subset(roads, RTTYP == "I")
roads <- st_transform(roads, crs = st_crs(counties.pre))

# Identify intersecting counties
intersection_list <- st_intersects(counties.pre, roads)
counties.pre$highway <- as.integer(lengths(intersection_list) > 0)

g <- counties.pre %>%
  st_drop_geometry() %>%
  dplyr::select(fips.pre, highway)

g$fips.pre <- as.numeric(g$fips.pre)

write_csv(g, output_file)
message("[OK] Saved highway data to: ", output_file)